% load default parameters for temperature simulation information
% all units are in m, this will break some previous code so fix as required
function [data,P] = load_sim_data(P0,r,z)

P.t_glass = 170;
P.r_total = 2000;
P.t_water = 2000;
P.w0 = 2;
P.r_fine1 = 20;
P.r_fine2 = 100;
P.r_fine3 = 300;
P.mesh1 = 1e-7;
P.mesh2 = 5e-6;
P.mesh3 = 1e-5;
P.mesh4 = 2e-5;
P.ins_bot = 1;  % whether bottom of glass plate is assumed to be insulating 
% or constant temperature to extrapolate behavior. Even though this is the
% source of cooling, since it's closest it might also dictate background
% temperature rises
% 0: fixed T, 1: insulating, 2: forced convection
P.u = 0.5;  % speed of convection
P.x = 0.01;
P.heatz = 3e-5;    % z position of water heating center

P0fields = fields(P0);
for i = 1:length(P0fields)
    eval(['P.',P0fields{i},'=P0.',P0fields{i},';'])
end

if nargin == 1
    zvals = P.heatz;
else
    zvals = [max(P.simz_list(P.simz_list<=z)),min(P.simz_list(P.simz_list>=z))];
end
temps = zeros(size(zvals));

for i = 1:length(zvals)
    fdn = 'COMSOL_simulation/';
    P.heatz = zvals(i);
    sim_file_name = [fdn,'vol_tg',num2str(P.t_glass),'_rt',num2str(P.r_total),'_tw',...
            num2str(P.t_water),'_w',num2str(P.w0),'_rs_',num2str(P.r_fine1),'_',num2str(P.r_fine2),...
            num2str(P.r_fine3),'_m',num2str(P.mesh1),'_',num2str(P.mesh2),'_',num2str(P.mesh3),...
            '_',num2str(P.mesh4),'_hz',num2str(P.heatz),'_u',num2str(P.u),'_x',num2str(P.x),'_',num2str(P.ins_bot),'.csv'];

    data = csvread(sim_file_name,9);
    T_base = 285.15;
    data(:,1:2) = data(:,1:2);
    data(:,3:6) = data(:,3:6) - T_base;

    if nargin > 1
        % interpolate to the specified r and z positions
        rs = [max(data(data(:,1)<=r,1)), min(data(data(:,1)>=r,1))];
        zs = [max(data(data(:,2)<=z,2)), min(data(data(:,2)>=z,2))];
        w1 = [1-(r-rs(1))/(rs(2)-rs(1)), (r-rs(1))/(rs(2)-rs(1))];
        w2 = [1-(z-zs(1))/(zs(2)-zs(1)), (z-zs(1))/(zs(2)-zs(1))];
        ind = zeros(2);
        temps(i) = 0;
        for ii = 1:2
            for jj = 1:2
                ind(ii,jj) = find(data(:,1)==rs(ii) & data(:,2)==zs(jj));
                temps(i) = temps(i) + data(ind(ii,jj),6)*w1(ii)*w2(jj);
            end
        end
    end
end

if nargin > 1
    data = ((z-zvals(1))*temps(2)+(zvals(2)-z)*temps(1))/(zvals(2)-zvals(1));
end